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Abstract 

We estimate generic statistical properties of a structural credit risk model by considering an ensemble of correlation 
matrices. This ensemble Is set up by Random Matrix Theory. We demonstrate analytically that the presence of correlations 
severely limits the effect of diversification In a credit portfolio if the correlations are not identically zero. The existence of 
correlations alters the tails of the loss distribution considerably, even If their average Is zero. Under the assumption of 
randomly fluctuating correlations, a lower bound for the estimation of the loss distribution Is provided. 
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Introduction 

The financial crisis of 2008-2009 clearly revealed that an 
improper estimation of credit risk can lead to dramatic effects on 
the world's economy. The vast underestimation of risks embedded 
in credits for the subprime housing markets induced a chain 
reaction that propagated into the worldwide economy. A better 
estimation of credit risk (see, eg, [1,2,3,4,5]) is therefore of vital 
interest. We can distinguish two fundamentally different ap- 
proaches to credit risk modeling (see, eg, [6]): the structural and 
the reduced-form approach. 

Structural models have a long history, going back to the work of 
Black and Scholes [7] and Merton [8] . The Merton model assumes a 
zero-coupon debt structure with a fixed time to maturity. The 
value of the company's assets is modeled by a stochastic process. 
The risk of default and the associated recovery rate, the residual 
payment in case of a loss, are directly determined by the 
company's asset value at maturity. 

Reduced-form models attempt to capture the dependence of default 
and recovery rates on macroeconomic risk factors. Both quantities 
are modeled as independent stochastic variables. Some well known 
reduced-form model approaches can be found in [9,10, 1 1,12,13]. 

First passage models were first introduced by Black and Cox [14] 
and they fall somewhat in between the two modeling approaches. 
Similar to Merton's model, the market value of a company is 
modeled by a stochastic process. However, in the first passage 
models a default occurs whenever this market value hits a certain 
threshold for the first time. The recovery rates are typically 
modeled independently, for example, by a reduced-form model, 
see [15,16], or are even assumed to be constant, see, eg, [6]. 
Recent approaches aim at improving first passage models by 
including the chance of full recovery, even if a company's market 
value is below the threshold, see [17], and estimating correlations 
between default probabihties of industry sectors, see [18]. 

Reduced-form and frrst-passage models are implemented in 
commercial software solutions, for example, CreditMetrks initially 
developed by JP Morgan [19], CreditPortfolioView by McKinsey & 



Company [20] or CreditRisk+ by Credit Suisse [21]. As there can 
be a strong connection between default risks and recovery rates, 
the chances of large losses are often underestimated in the 
reduced-form and first passage models, see [22,23]. The Merton 
model does not require this separation and is, for example, 
adopted by Moody's KMV. 

Structural models provide a "microscopical" tool to study credit 
risk as the defaults and recoveries are traced back to stochastic 
processes modeling the state of individual obligors. For a portfolio 
of credits, such as collateralized debt obligations (CDOs), 
correlations represent a key factor that influences its risk. The 
benefit of diversification, ie, the reduction of risk by increasing the 
portfolio size, is severely limited by the presence of even weak 
correlations. This has been demonstrated for the case of constant 
positive correlations, both in the first passage model with constant 
recovery [24,25] and in the Merton model [26,22]. The key 
problem in estimating the credit risk of a realistic portfolio is of 
course the huge number of parameters involved. This is precisely 
where approaches from statistical physics can be most helpful: the 
state of a system with many degrees of freedom is, under certain 
conditions, described by few macroscopic observables. In the 
thermodynamic equilibrium, these are energy, temperature, 
pressure, etc. Ergodicity holds, ie, time and ensemble average 
yield the same results. A somewhat similar situation exists for 
spectral statistics in quantum chaotic systems, see [27]. A moving 
average over one long spectrum equals an ensemble average over 
random matrices, if the number of levels is very large. Originally, 
random matrix theory was developed in the 1950s to describe the 
spectra of heavy nuclei, see [28] . Here we transfer this idea to large 
credit portfohos involving correlated assets. In the case of a great 
many contracts, we expect a self-averaging property which then 
should allow to average over an ensemble of random correlation 
matrices. We manage to carry out this approach largely 
analytically. We obtain estimates for the distribution of asset 
values and the portfolio loss distribution in which the complicated 
effects of all correlations are indeed reduced to a single parameter 
measuring the correlation strength. 
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A Structural Credit Risk Model 

Our model is based on Mertoii's original model, assuming a 
zero-coupon bond for the debt structure of the obligor. Our aim is 
to analytically describe the impact of correlations on the losses of a 
credit portfolio. Even though the Merton model makes many 
simplifying assumptions, it can provide more than just qualitative 
insights into credit risk. Indeed we demonstrated recently that 
empirical credit data are in accordance with analytical results 
derived from the Merton model [29]. 

The cash flow of the zero-coupon bond is limited to two dates: 
the date of issue t = 0 and maturity t=T. At the issue date the 
creditor lends a specified amount of money to the obligor. At 
maturity, the obligor has to repay the face value of the bond. The 
face value is the amount borrowed plus interest and risk premium. 
A default occurs if the asset value Vk of company k is below the 
face value at maturity time T. The size of the loss then depends 
on how far Vic is below the face value ft. We assume that the asset 
values in a portfolio of K companies foUow a geometric Brownian 
motion. An overview of the model's input parameters is given in 
Table 1. 



Equation (1) gives the pdf of asset values in the case of a 
correlated Brownian motion. However, we are not interested in 
the impact of a specific correlation matrix. Instead we want to 
estimate the general impact of correlations. To this end, we want 
to average over all possible correlation matrices and disclose the 
general statistical behavior. 

We use a random matrix approach to calculate the average 
distribution of asset values, ^*™*'(I^)), for random correlations 
where the average correlation level is zero. To achieve this we 
replace the covariance matrix S by 



l.w = SWW^S 



(3) 



where S = diag(f7i, . . . ,aii) contains the standard deviations and 
WeM^ ^ is a random matrix. The entries of W are independent 
and Gaussian distributed, 



Average distribution of asset values 

For the sake of simplicity, let us first consider the case of a 
Brownian motion for the asset values. Later on this can be easily 
mapped to the geometric Brownian motion by a simple 
substitution. For a Brownian motion, the probability density 
function (pdf) of the vector V of K asset values at maturity T is 
described by 




exp 



(V-nTyi.-\V-nT) 



(corr) 



exp 



N 



(4) 



with variance X/N. The resulting correlation matrix WW^ is 
Wishart-distributed [30] with average correlation zero. With the 
parameter A' we can control how strongly the entries of WW^ 
fluctuate. For N—kx>, we obtain the unit matrix for WW\ ie, the 
uncorrelated case. For N>K, we obtain an invertible covariance 
matrix with random entries. The case N <K is disregarded as the 
resulting matrix is not invertible which is usually required for 
applications in risk management. When mserting this ansatz into 
Eq. (2), we obtain 



Here, S is the covariance matrix and fi is the drift vector. For 
later convenience we can express this as a Fourier transform. 



</;("^)(F)> = 



p'^''°"\w)p^'"^'i{V,SWW^S)d[W] (5) 



1 



(271)- 
exp 



exp ( — io3\V — 12 T)) 



-oj'Sco )d[co 



(2) 



i2nf J 



exp(-/a;tK)— ^d[w] (6) 

^det(NI+TSo}w^S) 



Table 1. Input of the structural credit risk model. 





Variable 


Description 


Unit 


K 


Number of contracts 




T 


Time to maturity 


[year] 


"k 


Volatility of asset k 


[year]-^'^ 


V-k 


Drift of asset k 


[year]-' 


N 


Parameter to control correlations, 




/y—^co; uncorrelated limit 




Initial value of asset k 


[currency] 


h 


Face value of contract k 


[currency] 



doi:l 0.1 371 /journa!.pone.0098030.t001 
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Figure 1. Illustration of the average asset value distribution (jf'^''\py) for T = 1, K = 50 and different values for N. Solid, 
dashed{dotted, dashed and dotted lines correspond to N = K, 2K, 5/< and 30K, respectively. 
doi:1 0.1 371 /journal.pone.0098030.g001 



where / denotes the unit matrix. A detailed derivation is given in 
appendix SI. Here we choose ^ = 0. We will reintroduce the drift 
later on, when we make the substitution for the geometric 
Brownian motion. The determinant can be written as 



(7) 



K 



2nf r(N/2) l^A = i oi, 

N-K 



N 



K. N-K I 



fll) 



because the matrix S(oo3^ S has rank one. Hence, we arrive at 



<P""^'(K)> = 



(271)- 



nf\ 



exp( — /ffl'K) 



(\+{TlN)oi^SScofl^ 



(8) 



for the average distribution of if assuming a randomly 

distributed correlation matrix and an average correlation level of 
zero. We stated earlier that we include N in the distribution of the 
random matrices W in order to render the variance of the average 
asset value distribution A'-independent. The variances only 
depend on T and cr^-, as discussed in appendix S3. The parameter 
N is only used to control the correlations. In hyperspherical 
coordinates, Equation 1 1 depends only on the hyperradius 



This integral can be calculated by using the Gamma function 
(see [31]) in the form 



z' 'exp( — az)dz, x>0, a>0 



(9) 



We identify a"'^ with {\ +{T / N)o3' SSwY^ and obtain 



1 



1 



■exp(- 



exp 



N V? 



ATz- 



vE^ (10) 



as worked out in appendix S2. This integral is a representation of 
the Bessel function of the second kind K, of the order {K — N)/2, 
see [32]. Thus, we obtain 



This leads to the expression 
<P""^*(P)> 



-K 



1 

2 2 N+K-l 
P 2 



2tiT T(N/2) 

N-K 



(12) 



(13) 



for the hyperradial density function, cf appendix S3. We illustrate 
this density function in Figure 1 for .^T = 50 and different values of 
N . We note that the tail-behavior for large p is exponential. 

We obtain the average asset value distribution in case of a 
geometric Brownian motion by a simple substitution Fi;-> Vk, 
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-K 1 N 

N 2'~2 



2kT r(N/2) \k=i ffkVk 



K 

n 



K, N-K 
2 I 




(14) 



with 



(15) 



Here, the parameter refers to the standard deviation of the 
underlying Brownian motion, ie, the volatility of asset returns. The 
resulting asset values thus have the variance 



a^=exp(2/.^r)(exp(4r)-l)Ki 



(16) 



where Vkfi are the starting asset values at ? = 0. Figure 2 shows the 
distribution of asset values based on a geometric Brownian motion, 
as given in Eq. (14). The findings are similar to the case of the 
Brownian motion. While we obtain a narrow but heavy-tailed 
distribution for N = K, the distribution slowly approaches an 
uncorrelated bivariate log-normal distribution with increasing 
values of A'^. 



Lk^ 



Vk<Fk 



0 , else 



(default) 
(no default) 



(17) 



We observe that the asset values have to be positive in Eq. (17). 
Therefore we assume in all further considerations that the 
underlying asset value process is given by a geometric Brownian 
motion. 

When calculating the overall loss of a portfolio, we have to 
weight each loss by its face value in relation to the sum of all 
portfolio face values, 



L = '^fkLk , fk 



Fk 



(18) 



We integrate over the pdf of asset values and filter for those that 
lead to a given total loss L. By the above stated definitions, we can 
define a filter for the total loss at maturity time T. In the next step 
we express the filter using a Fourier transformation. Eventually, we 
separate those terms that correspond to a default and those that 
describe the asset values above the face value F^. 



„(loss) 



(L)- 



00 



iV)S\L 



(19) 



Loss distribution 

We now turn to the calculation of the loss distribution. A default 
occurs if the asset value F<. at maturity T is lower than the face 
value ffc. The size of the loss is given by the difference of F/t and 
Vk- Even if a loss occurs, the creditor might not lose aU money that 
he lent, because the obligor is still able to pay back the amount Vk. 
In order to compare losses in a portfolio of credits, we have to 
normahze them by the corresponding face value. We define the 
normalized loss L/^ of the k-ih asset as 



1 

271 



dv exp{—ivL) 



d[ Flijf^^'C V) — dv exp ( - ivL + iv J^fkLk ] (20) 



d[F]exp|^/vg/,L,jy-^'(K) (21) 




Figure 2. Illustration of the average asset value distribution </)''"^'(F))with a geometric Brownian motion for T = ^, K = 2, 1^0 = 
100, /(= 0.05 and different values for /V. Both distributions have the identical standard deviation ax 16 ((r = 0.15). For W = 2, we obtain a heavy- 
tailed distribution while the uncorrelated limit is reached for N = 100. 
doi:1 0.1 371/journal.pone.0098030.g002 



PLOS ONE I www.plosone.org 



4 



May 2014 I Volume 9 | Issue 5 | e98030 



A Random Matrix Approach to Credit Risk 



with 



K 
k=\ 



1 

In 



dFfcexp ivfki 1 



dv exp( — ivL) 



+ 



ymv)(^ (22) 



r(v,z)= n 

k = l 



d Vk exp ( ivfk — — I + 



dVk 



2o't VkV nzT 



. /N_l_NilniVk/Vk,o)-(^k-4/2)Tf\ (^5) 



Here, the expression in the square brackets acts as an operator, 
because does not necessarily factorize. We will use this 

ansatz to calculate the average loss distribution in the next section. 
However, Eq. (22) can be used to calculate the loss distribution if 
the actual asset value distribution is known, ie, the statistical 
dependence and the underlying process are estimated. To prepare 
for this, it is handy to write Eq. (22) as a combinatorial sum. 



^(ioss)(^-)_ 1 [ dvexp(-ivZ,) (23) 
2n J 



We carry out a second order approximation of this expression in 
appendix S4 and arrive at 



1 

'2nr(N/2) J 



VW2(Z) 



exp 



dz exp( - 



(L-mi(z)y 
2m2{z) 



(26) 



with 



K (fc) 

E E 

k=l j=l 



n 



j dV,sxp(ivf,^) 



/ePenn(/,A:,^ 

00 

n s^^i 

qe{l...K} ^? 



wi (z) = '^fkmi,k(z) (27) 



^2(^) = J2f^(m,k(^) - miM) (28) 



where Ve,Tm(j,k,K) is the y'-th permutation of fc elements of the set 
{ 1 ... AT} . For example, if K=3 and k = 2, we obtain, 
Perm(l,2,3) = {l,2}, Perm(2,2,3) = {2,3} and Perm(3,2,3) = 
{1,3}. However, Eq. (24) might need to be estimated numerically, 
depending on the complexity of the asset value distribution 
p^'^'"\V). In the section Homogeneous portfolios, we wiU simplify 
this combinatorial sum for a homogeneous portfoUo and the 
average asset value distribution <p<"'''(I^)>. 

Average loss distribution 

Now we have developed all necessary tools to model the average 
distribution of losses, under the assumption of random correlations 
and an average correlation level of zero. We start by inserting the 
average asset value distribution in a component-wise notation (cf. 
appendix S2) into the loss distribution (22), 

00 

I dvexp( — (vZ,)j-(v,z) 

— 00 



and 

m i,{z)= ^ [ 1- ( Pk-Vk V 
2ak\/nzT \ Vk\ Fk J 
0 

xexp(-^(^"^^^/W..-c^l/2)r)^j,^^ (29) 

However, the convergence radius of the power series expansion 
involved in this approximation is one. Although we consider large 
portfolios K,ie,fk is small, v runs from — 00 to + co. This second- 
order approximation might describe the default terms adequately. 
However, the non-default terms, corresponding to a delta peak at 
L = 0 require v to run from — 00 to + 00 . Thus, the non-default 
terms cannot be approximated using this second-order approxi- 
mation. To circumvent this problem we develop an improved 
approximation in the next section. 

Due to the complexity of wi(z) and »!2(z), the z integral needs 
to be evaluated numerically. We present this for the example of a 
homogeneous portfolio. 
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Homogeneous portfolios 

In case of a homogeneous portfolio, in which all credits have the 
same face value Fk = F and the same variance ff^ = ff^ and initial 
value Fi,o = Vo, the weights can be simplified to 



fk 



K 



(30) 



As m\ji{z) and wi,t(z) become identical for every k, we denote 
them by fni(z) and m\{z) leading to 



1 1 ( ^{\HFIVo)-(ji-o'/2)T) \ 

2 2 \ 2aV?T ) ^ ' 



and the default term (r*'^'(v,z)y where 

F 

r(°)(v,z)=|dKexpg 



F-V 



Wi(z) = wi(z) 



1 2 
W2(z)= ^(W2(z)-Wl(z) ) 



(31) 



(32) 



W;(z) = 



//V f 1 fF-V 



2(j^/nzf} V\ F 

0 



(33) 



Here F is a scalar and we only have to calculate a single integral 
over V. After inserting this into Eq. (26), we can calculate the loss 
distribution for a homogeneous portfolio in the second order 
approximation. 

Improved approximation for a homogeneous portfolio 

The second order approach can be improved by approximating 
the indi\'idual terms of the loss distribution instead of approxi- 
mating the expression as a whole, similar as discussed in [26]. In 
case of a homogeneous portfolio the combinatorial sum in Eq. (24) 
reduces to 



GO -t- 00 



7=0 ^ J 

with the non-default term (r'^''') where 



(34) 



„(ND). 



dV- 



X exp 

2(7 VynzT 



Ni\niV/Vo)-iii-ay2)Tf ^ 
AzTa^ , 



(37) 



In the homogeneous case, the integration variable F is a scalar. 
The approximation follows the same principles as in the previous 
section, resulting in 



-t-co 

I dvexp(-(vZ.)(r<°)(v,z)y 



-hCC 

= 1 



dvexp( iv[ ^mi(z)-Lj - ^ (m2{z)-mi{zf) ) (38) 



y(w2(z)-mi(z) 



-^^^^expf-#^=/!!^| (39) 



2/'(w2(z)-wi(z) 



In this approximation, the non-default terms given by Eq. (36) 
can be calculated exactly. They correspond to a delta peak at 
L = Q. Another advantage over the approach presented in Eq. (26) 
is that the approximation is performed for each number of defaults 
j separately and weighted by j/K accordingly. Here, the omitted 
third term is of the order j/K^ and thereby much smaller than the 
third term of the simple second order approximation (33), which 
would be of the order \jK} . Thus, when approximating each term 
in the combinatorial sum separately, we obtain an improved result. 
Insertion into (34) leads to 



1 

2nT(N/2) 



7=0 ^ i 



dz z2 exp(— z) 



2aVy/nzT 



\j(m2{z)-m\{zf) \ 2j(rn2(z)-m\{zf] 



xexp,- ^(^-^^/^;^r-W ) (35) 



x(i + iErf 



MFlVo)-(ii-a^l2)T 



2(j\/zT 



(40) 
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Figure 3. The loss distribution iot K = 10, >7 = OAS,// = 0.05, T 
= ^, Vo = ^00, F = 75 and different amounts of randomness in 
the correlation matrix, N = K (solid black), N = 2K (dashed 
blue), N = 10/r (dotted red), N = 30/(^ (dot-dashed green). 

doi:1 0.1 371 /journal.pone.0098030.g003 

which is the final result. 

Results 

We now apply the analytically developed model to a specific 
example. To analyze the impact of correlations, we calculate the 
loss distribution for different homogeneous portfolios with sizes 
^=10, ii: = 50 and K=\00 with the parameters Fo = 100, 
;i = 0.05, f7 = 0.15, F = 75 and T=l. As stated in the previous 
section, we can control the amount of correlation in our model 
with the parameter N. Since we only consider correlation matrices 
with full rank, we obtain the strongest correlations if we choose 
N = K. For N—*<X), the correlation matrix becomes the unit 
matrix. Thus, this represents the transition to a system without 
correlations. As we have to evaluate the loss distributions 
numerically, the limit N^oo has to be properly interpreted. We 
need to identify a value for which this convergence is valid in good 
approximation. Figure 3 illustrates the loss distribution for K=W 
and different values of N. Our study indicates that a value of 
N = 30K is a good choice for approximating the uncorrelated case 
and is stiU numerically feasible. The results are presented in Fig. 4. 
For all portfolio sizes, K=10, K = 50 and .ST =100, we obtain 
heavier tails of the loss distribution of the correlated portfolio 
compared to the uncorrelated case. Even the simple approxima- 
tion, represented by the dashed blue curve, exhibits these heavy 
tails. With the inserted logarithmic plots, we can identify a nearly 
power-law decay of the loss distribution for the correlated case. 

The distributions become narrower for larger values of K. 
However, the tails of the correlated case remain heavier than those 
of the uncorrelated case. While both approximations yield similar 
results for K= 10, their difference becomes larger with K. As both 
approximations have to be performed numerically, the improved 
approximation is always favored. However, the tail behavior 
remains the same, even for the simple approximation, as indicated 
by the logarithmic scaled inserts in Fig. 4. This is a strong 
indication that the tails of the loss distribution are vastly 
underestimated if correlations are not taken into account. 

Due to the approximation, the normalization of the loss 
distribution is not exact. Especially the normalization of the 
simple approximation is problematic for large values of K. The 
normalization might also be used as an indication for the quality of 
the approximation. The improved approximation exhibits a delta 
peak at L = 0, as the non-default terms can be calculated exactly. 
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(a) K=10 




(b) K=50 




(c) K=100 

Figure 4. The loss distribution of a homogeneous portfolio 
with a- = 0.15, fi = 0.05, T = 1, Uo = loo, F = 75 and different 
values of K. The blue dashed line represents the simple approxima- 
tion; the solid black line represents the improved approximation. Both 
have been calculated for the strongest random correlations, N = K. The 
uncorrelated case is given by the red dotted line, calculated with the 
improved approximation with N = 30/C. 
doi:10.1371/journal.pone.0098030.g004 

However, the interval [0; 0.0002[ was not evaluated due to 
numerical feasibility. 

In our example, we do not vary the maturity time T, ie, we 
choose 7" = 1 . One can increase T to estimate the evolution of the 
loss distribution. However, this evolution depends strongly on the 
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drifts and standard deviations at - Depending on their value, the 
exposure to default risk can either increase or decrease. 

Discussion 

To assess the risk of a credit portfolio, it is crucial to take 
correlations between obligors into account. We consider the 
Merton model, in which defaults and recoveries are determined by 
the underlying asset processes. The correlation matrix of the asset 
returns has to be estimated from historical time series. This is not 
always easy, because the correlations change in time, ie, they are 
non-stationary. Since only time series up to a certain length can be 
used, the correlation coefficients contain a specific t^-pe of 
randomness, see [33,34]. Several methods have been put forward 
to estimate and to reduce this "noise". Thus, we assume that such 
a noise reduction has been done. The corresponding "true" 
correlation coefficients and matrices are the proper input for the 
structural credit risk model of the Merton t)'pe that we consider. 
We discussed this issue of noise reduction to emphasize that the 
random matrix approach in that context focuses on the spectral 
statistics of correlation or covariance matrices, see 
[35,36,37,38,39]. It is based on a very different motivation as 
compared to the present application. 

Searching for generic properties, we devised the present random 
matrix approach. Instead of calculating the portfolio loss 
distribution for a specific correlation matrix, we average over an 
ensemble of random correlation matrices. Our approach transfers 
concepts of statistical physics. In quantum chaos, the average over 
an individual, long spectrum equals the average over an ensemble 
of random matrices, if the level number is very high. We expect 
that a similar self-averaging property also holds here. This line of 
reasoning is supported by the following consideration: The 
correlation coefficients are varying functions in time, because the 
business relations of the companies change. This impUes that a 
correlation matrix over a somewhat longer period in time is a 
varying quantity, ie, it corresponds to some kind of ensemble. 

In our model the average correlation level is zero and we 
assume that there is no branch structure in the correlations. The 
fluctuation strength of individual correlations is controlled by a 
single parameter. This ansatz allowed us to estimate generic 
statistical properties of the Merton model. Some features are not 
taken into account which are present in empirical data, such as 
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